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A method for enhancing the stability and robustness of explicit schemes in compu- 
tational fluid dynamics is presented. The method is based in reformulating explicit 
schemes in matrix form, which can then be modified gradually into semi or strongly- 
implicit schemes. 

Prom the point of view of matrix-algebra, explicit numerical methods are special 
cases in which the global matrix of coefficients is reduced to the identity matrix /. 
This extreme simplification leads to severe limitation of their stability range, hence 
■ of their robustness. 

In this paper it is shown that a condition, which is similar to the Courant- 
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Friedrich-Levy (CFL) condition can be obtained from the stability requirement of 
inversion of the coefficient matrix. This condition is shown to be relax-able, and that 
a class of methods that range from explicit to strongly implicit methods can be con- 
structed, whose degree of implicitness depends on the number of coefficients used 
in constructing the corresponding coefficient-matrices. Special attention is given 
^ ■ to a simple and tractable semi-explicit method, which is obtained by modifying 

the coefficient matrix from the identity matrix I into a diagonal- matrix D. This 
method is shown to be stable, robust and it can be applied to search for station- 
ary solutions using large CFL-numbers, though it converges slower than its implicit 
counterpart. Moreover, the method can be applied to follow the evolution of strongly 
time-dependent flows, though it is not as efficient as normal explicit methods. 

In addition, we find that the residual smoothing method accelerates convergence 
toward steady state solutions considerably and improves the efficiency of the solution 
procedure. 
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1 Introduction 



In the last two decades, tremendous progress has been made in both compu- 
tational fluid dynamics (CFD) algorithms and the computer hardware tech- 
nologies. The computing speed and memory capacity of computers have in- 
creased exponentially during this period. Similarly is astrophysical fluid dy- 
namics (AFD), which is considered nowadays as a rapidly growing research 
field, in which modern numerical methods are extensively used to model the 
evolution of rather complicated flows. Unlike CFD, in which implicit meth- 
ods are frequently used, the majority of the methods used in AFD are ex- 
plicit. Several of them became very popular, e.g., ZEUS (Stone & Norman 
1992), NIRVANA+ (Ziegler 1998), FLASH (Fryxell et al. 2000), VAC (Toth 
et al.1998), THARM (Gammie et al. 2003). The popularity of explicit meth- 
ods arises from their being easy to construct, vectorizable, parallelizable and 
even more efficient as long as the dynamical evolution of compressible flows 
is concerned. Specifically, for modeling the dynamical evolution of HD-flows 
in two and three dimensions explicit methods are highly superior to-date. For 
modelling relativistic flows, Koide and collaborators (e.g., Koide et al. 2000, 
2002, Meier et al. 2001) and Komissarov (2001) have developed pioneering 
general relativistic MHD solvers. A rather complete review of numerical ap- 
proaches to relativistic fluid dynamics is given by Marti & Miiller (1999) and 
Font (2000). A ZEUS-like scheme for general relativistic MHD has also been 
developed and is described in De Villiers & Hawley (2003). 
However, these methods are numerically stable as far as the Courant-Friedrich- 
Levy number is smaller than unity. The corresponding time step size decreases 
dramatically with the incorporation of real astrophysical effects. Specifically, 
they may even stagnate if self-gravity, radiative and chemical effects are in- 
cluded. Moreover, explicit methods break down if the flow is weakly or strongly 
incompressible, and if the domain of calculations is subdivided into a strongly 
stretched mesh. In an attempt to enhance their robustness, several alterna- 
tives have been suggested, such as semi-explicit, semi- implicit or even implicit- 
explicit methods (Kley 1989, Toth et al. 1998 ). Nevertheless, their rather 
limited range of applications has lead to the fact that most of the interesting 
astrophysical problems remained, indeed, not really solved. A simple example 
is the evolution of a steady turbulent accretion disk. It was found by Balbus 
& Hawley (1991) that weak magnetic fields in accretion disks are amplified, 
generate turbulence, which in turn redistribute the angular momentum in the 
disk. However, whether this instability leads to the long-sought global steady 
accretion rate, or is it just a transient phenomenon in which the generation of 
turbulence is subsequently suppressed by dynamo action are not at all clear. 
Other notable phenomena are the formation and acceleration of the observed 
superluminal jets in quasars and in microquasars, the origin of the quasi- 
periodic oscillation in low mass X-ray binaries or the progenitors of gamma 
ray burst are still spectacular. 
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Fig. 1. Schematic description of several coefficient matrices used by different numer- 
ical solution procedures in two dimensions. Most commonly used matrices are Fl, 
F2, F3, F4 and F5, which correspond to tri-diagonal block, scalar tri-diagonal, block 
diagonal, scalar diagonal and identity matrices, respectively. Solution methods that 
rely on the inversion of Fl are classified as implicit, whereas those relying on F5 are 
explicit. F2, F3 and F3 correspond to intermediate methods, i.e., explicit-implicit 
methods. 
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Fig. 2. A schematic description of the time step size and the computational costs 
versus the band width M of the coefficient matrix. N is the number of unknowns. 
Explicit methods correspond to M = 1 and large 1/St. They require minimum 
computational costs (CC). Large time steps (i.e., small l/5t) can be achieved using 
strongly implicit methods. These methods generally rely on the inversion of matrices 
with large band width, hence computationally expensive, and, in most cases, are 
inefficient. 

Explicit methods rely on time-extrapolation procedures for advancing the so- 
lution in time. However, in order to provide physically consistent solutions, it 
is necessary that these procedures are numerically stable. The usual approach 
for examining the stability of numerical methods is to perform the so called 
von Neumann analysis (see Hirsch 1988 for further details). For example, using 
an explicit procedure to solve the simple advection-diffusion equation: 



and applying the von Neumann stability analysis, it is necessary, but not suf- 
ficient, that the CFL- number is less than unity. This is equivalent to require 
that the time step size fulfills the inequality: 5t < min( WTTTi i ^7~)> wnere 

T, is, f, p denote an arbitrary diffusible variable, viscosity coefficient, internal 
or external forces and density of matter, respectively. The sub-scripts t, x, xx 
represent first and second order derivatives of T with respect to time and to 
the x-coordinate. The minimum function in the last inequality runs over each 
grid point of the domain of calculation. For additional details see Sec. 9.4 in 
Fletcher (1991). The force / may corresponds to internal forces, such as ther- 
mal, radiative or magnetic pressure. In regions where the density is low (e.g., 
corona around compact objects such as black holes or even above accretion 
disks), the derivative |f^| 1//2 corresponds to the sound speed, speed of light 
or to the propagational speed of magnetic Alfven waves, all of which can be 
extremely large compared to the actual velocity of the flow. Consequently, the 
CFL-condition limits the range of application and severely affects the robust- 
ness of explicit methods. In particular, equations corresponding to physical 
processes occurring on much shorter time scales than the hydro-time scale 



T t + uT x = z/T xx + /, 
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(e.g., radiation, self-gravitation and chemical reactions) cannot be followed 
explicitly. Furthermore, these methods are not suited for searching solutions 
that correspond to evolutionary phases occurring on time scales much longer 
than the hydro-time scale. Using high performance computers to perform a 
large number of explicit time steps may lead to accumulation of round-off 
errors that can easily distort the propagation of information from the bound- 
aries and cause divergence of the solution procedure, especially if Neumann 
type conditions are imposed at the boundaries. 

In contrast to explicit methods, implicit methods are based on solving a ma- 
trix equation of the form A5q = d, where A is the coefficient matrix resulting 
from the linearization of the system of equations to be solved, d is the right 
hand side vector of known quantities, and Sq is the solution vector sought. 
These methods have two major drawbacks. First, constructing the matrix A 
is difficult, time consuming, and may considerably influence the stability and 
robustness of the method. Second, the inversion procedure must be stable and 
extremely efficient. In general, conservative discretization of the MHD equa- 
tions give rise to sparse matrices, or even to narrow band matrices. Therefore, 
any efficient matrix inversion procedure must take the advantage of A being 
sparse (see Fig. 1). Inverting A directly by using Gaussian elimination requires 
N 3 algebraic operations, where N is the number of unknowns (Fig. 2). If the 
flow is multi-dimensional and a high spatial resolution is required, the num- 
ber of operations can be prohibitive even on modern supercomputers. Krylov 
sub- iterative methods (KSIMs), on the other hand, are most suited for sparse 
matrices and avoid the fill-in procedure. In the latter case, A is not directly 
involved in the process, but rather its multiplication with a vector. The conver- 
gence rate of KSIMs has been found to depend strongly on the proper choice of 
the pre-conditioner. For advection-dominated flows, incomplete factorization 
such as ILU, IC and LQ, approximate factorization, ADI, line Gauss-Seidel are 
only a small sub-set of possible sequential pre-conditioners (see Saad & van 
der Vorst 2000 and the references therein). Another powerful way of accelerat- 
ing relaxation techniques is to use the multi-grid method as direct solver or as 
pre-conditioner (see Brandt 2001; Trottenberg et al. 2001). For parallel com- 
putations, Red-Black ordering in combination with GRMES and Bi-CGSTAB 
as well as domain decomposition are among the popular pre-conditioners (see 
Dongarra et al. 1998 for further discussion). 

In this paper we present a strategy for enhancing the stability and robustness 
of explicit methods. 

In Sec. 2 of this paper we describe and apply the new strategy to scalar equa- 
tions in one-dimension, and generalize it to system of equations in Sec. 3. 
Here it is shown that a class of semi-explicit numerical methods can be eas- 
ily constructed through reformulating explicit methods in matrix form, which 
thereafter can be modified. Spacial attention is given to a simple semi-explicit 
method which can be applied for searching stationary solutions. Although The 
method is easy to program and stable even when using CFL > 1, it converges 
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relatively slowly compared to implicit methods, and therefore further addi- 
tional acceleration techniques of convergence are still to be found. The results 
of various test calculations are presented in Sec. 5, followed by the summary 
and conclusions in Sec. 6. 



2 Stability of matrix inversion and the CFL connection 



In fluid flows, the equation of motion which describes hydrodynamically the 
time-evolution of a quantity q in conservation form reads: 

ff + L# = /, (2) 



where V and / are spatially varying velocity field and external forces, respec- 
tively. L represents a first and/or second order linear differential operator that 
describe the advection and diffusion of q. In the finite space 3ft, we may replace 
the time derivative of q by: 

5q _ g D+1 - g D 

St~ St ' [) 



where q n and q n+1 correspond to the actual value of q at the old and new time, 
levels, respectively. 

An explicit formulation of Eq.3 reads: 

| = [-LqV + JT, (4) 



whereas the corresponding implicit form is: 

Sq 
St 



[-LqV + fr 1 . (5) 



Combining these two approaches together, we obtain: 

= ~Jt + e[ ~ Lq ^ + /r+1 + (1 " *)H^ + ^° = RHS > ^ 



where 6(0 < 9 < 1) is a switch on/off parameter. 

The spatial operator H = [—LqV + /] should be computed at each grid point, 
using a conservative discretization. In this case, the coefficients of the Jacobian 
matrix can be constructed by computing <9RHS/<9g n+1 which we decompose 
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in A = —[dH/dq} n+1 and I /St. Thus, Eq. 6 can be replaced by the matrix- 
equation: 

+ 6A)5q = RHS. (7) 

ot 



As a starting condition in each time step, it is suggested to set q n+l = q n . 
Thus, the scheme degenerates into explicit scheme for 6 = and into iterative 
implicit for 6 = 1. 

Eq. 7. is a combination of the Newton and the Crank-Nicolson methods 
(Fletcher 1991; p. 165 and 304). It requires that in each time step the RHS 
should vanish completely. If steady state or quasi- stationary solutions are 
sought, then both the RHS as well as (q n+1 — q n )/5t must vanish simulta- 
neously. In the latter case, care should be taken to assure that 6 — » 1 as 
6t — > oo. 

We note that since the LHS and the RHS of Eq. 7 are generally different, a 
loop of iteration within each time step must be constructed in order to obtain 
a reasonable value for 2 q n+1 . 
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Fig. 3. The neighboring block matrices in the x and y-directions resulting from 5-star 
staggered grid discretization. Entries marked with 'X' denote the elements usually 
used in the implicit operator splitting approach (Hujeirat & Rannacher 2001), '*' 
and '+' are coefficients corresponding to the the source terms. The semi-explicit 
method for a scalar equation relies on inverting the diagonal matrix whose entries 
are marked with X surrounded by squares. The generalization of the semi-explicit 
method to the multi-dimensional HD-equations requires inverting the block diagonal 
matrix -D mo d- 

The matrix A contains coefficients such as V/Ax, rj/Ax 2 that correspond to 
advection and diffusion terms, respectively, as well as additional coefficients 
that correspond to source terms. 

2 Such deviations may evolve if the LHS of Eq. 7 is not precisely the real Jacobian 
of Eq. 6, and if the equations to be solved are partially non-linear . 
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Noting that Eq. 2 describes the time evolution of a real physical variable q, and 
applying a conservative discretization to the advection operator L, e.g., first 
order up-winding, it is easy to verify that the resulting matrix A — I + 5t0A 
is positive definite (nonsingular) and, diagonally dominant. Therefore, the 
matrix A can be stably inverted (Hackbusch 1994). 

In this case the relevant question is: how far can we simplify A, and still able 
to provide corrections Sq, appropriate for convergence ? 
Let us replace the real matrix coefficient A by a matrix A, which is simple 
and easy to invert. Thus, instead of solving the "difficult" matrix equation 
A5q = 5tRHS, we solve A5q = 5tRHS. However, in order to assure that 
we are dealing with the original problem still, the matrices A and A must 
share the essential spectral properties, i.e, the "spectral equivalence" (Golub 
& Loan 1989; Hackbusch 1994, p. 217). For example in explicit methods A = I. 
In order that the matrices / and A be spectrally equivalent, the sum of the 
absolute values of the elements in each row of St 9 A must be smaller than the 
corresponding diagonal element of I. This implies that 2 |V|/Ax < 1/St. The 
latter condition is relatively strong, and a weaker condition can be obtained if 
the flux difference, rather than the flux itself, is considered. In this case, and 
in the absence of diffusion and external forces, the following inequality holds: 

\AVq\/Ax < \q\/5t, (8) 

where q may acquire negative and positive values. However, since Vq is a con- 
servative quantity, we may obtain an upper limit for the flux-difference: 

\AVq\ < max{(Vq) in , (Vq) ont } = \V\\q\, (9) 

where the sub-scripts "in, out" denote the locations of the in- and out-flows 
through the surfaces of an arbitrary finite volume cell. In writing the last 
equality, we have omitted these subscripts for simplicity. 
Thus, in terms of Eq. 7, neglecting OA is equivalent to require: 

\AVq\/Ax < |V||g|/Aa; < \q\/St. (10) 

This is equivalent to the outcome of the von Neumann stability analysis (see 
Richtmyer & Morton 1967), which yielded the well-known condition CFL = 
\V\5t/Ax < 1. 

It should be noted, however, that the classical derivation of the CFL con- 
ditions relies on the assumption that the velocity V in Eq. 2 is constant. 
Practically, the HD-equations are non-linear, and using CFL = 1 does not 
prevent the exponential growth of the instability. Indeed, most of the explicit 
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Fig. 4. The evolution of the residual and the CFL-number versus number of iter- 
ation for the heat diffusion equation (Eq. 21). In (a) the time step is constructed 
directly from the residual. Here the diagonal elements, or equivalently the time step 
size, oscillates around a canonical value that corresponds to CFL = 1. The matrix 
inversion is highly unstable and fails to provide a perturbation-free steady solution. 
In (b) the time step is set to increase gradually from a small value up to a value 
that corresponds to CFL=1.75. The residual in this case grows exponentially, which 
implies that the matrix coefficient cannot be stable-inverted. In (c) the diagonal 
elements, i.e., I /St, are taken to correspond to CFL=0.975. In this case the matrix 
has a stable inversion, hence converges, though extremely slow. 

methods used in astrophysical fluid dynamics provide stable solutions, if only 
the CFL-number is strictly less than unity, and in most CFL < 0.6 

is required. The matrix A can be decomposed as follows: A = D + L + U, 
where D is a matrix that consists of the diagonal elements of A. L and U con- 
tain respectively the sub- and super-diagonal entries of A. Noting that a con- 
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servative discretization of the advection-diffusion hydrodynamical equations 
(Navier-Stokes equations) gives rise to a D which contains positive values, we 
may reconstruct a modified diagonal matrix D mod — I/St + 9D. In this case, 
Eq. 7 is equivalent to: 

[D mod + e(L + U)]5q = RHS. (11) 

A slightly modified semi-explicit form can be obtained by neglecting the entries 
of the matrix 6(L + U). In this case, a necessary condition for the iteration 
procedure to converge is that the absolute value of the sum of elements in 
each raw of 9(L + U) must be much smaller than the corresponding diagonal 
element of D mod in the same raw. In terms of Equation 9, the method is said 
to converge if the entries in each row of -D mo d fulfill the following condition: 

l/5t + 0(\V\/Ax + r)/Ax 2 + g) > \\A - D^, (12) 

where \\A — D||oo denotes the oo— norm of A — D, i.e., the maximum row 
sum of the modulus of the elements of A — D. This condition can be fulfilled, 
however, if the flow is smooth, viscous, and if appropriate boundary conditions 
are imposed 3 . Consequently, the inversion process of D^^Sq = RHS should 
proceed stable even when large CFL- numbers are used (see Fig. 6 and 7). 

As a numerical test, we have applied this approach to the one-dimensional 
diffusion equation in spherical geometry (see Eq. 21, see Sec. 5.1, for a detailed 
description of the physical problem). 

In this particular model, the matrix A of Eq. 7 is removed, and the equation 
to be solved reads: 

(^)5q = RHS. (13) 

As a first step, the equation is solved using a direct band matrix solver, where 
the band width is set to 1. Thus, the solver in this case is exact up to the 
machine accuracy. The time step size is set to be determined from the resid- 
ual, without knowing a-priori about the CFL associated-problems. The cor- 
responding results are displayed in Fig. 4. Indeed, we see that the optimal 
values of the diagonal elements that limit the exponential growth of the resid- 
ual oscillate around a canonical value which corresponds to CFL=1 (a/Fig 5). 
Decreasing the diagonal values artificially, so to correspond to CFL number 
that is slightly larger than one, gives rise to an exponential growth of the resid- 
ual (b/ Fig. 4). On the other hand, setting the diagonal values to correspond 

3 Note that diffusion pronounces the inequality in Eq. 12, which gives rise to larger 
CFL- numbers. 
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Fig. 5. In the finite volume discretization, the temperature Tj is defined at the cell 
center r™, whose boundaries are rj rj+±. 

to CFL-numbers that are slightly smaller than unity, we find that the matrix 
can be stable-inverted, and the resulting solutions converge to the stationary 
solution, though extremely slowly (c /Fig. 4). 

To examine the connection between the stability of the matrix inversion and 
the CFL condition applied to the heat diffusion equation in dt, we may re-write 
this equation using finite volume descritization (Fig. 5), and use the defect- 
correction formulation of Eq. 7. Thus, at an arbitrary finite volume cell j, the 
equation reads: 

[S^STj^ + [1/St + DjjSTj + [S J+1 ]ST J+1 = RHS 3 , (14) 



where 5 J -_ 1 (= -[u/volj] x [rJ/Ar™]), S j+1 (= -[v/volj] x [r] +1 /Ar™ and 
Dj(= —S_j-i — Sj + i) correspond to the sub-diagonal, super-diagonal and to the 
diagonal entries of the coefficient matrix. Here ArJ 1 = {rj l _ l —rj l ), volj = (r| — 

r| +1 )/3 and RH Sj = -(^^) + (^){r?(^) - r? +1 (^ + \ 

J J 3 + 1 



Let us consider the following three cases: 

(1) The matrix corresponding to Eq. 14 is diagonally dominant for every time 
step size St, hence it can be stable-inverted. 

(2) If we neglect the sub- and super-diagonal entries and preserve Dj, then 
the diagonal dominance of the matrix is not affected, maintaining thereby 
its stability for any time step size St. 

(3) If we decide, however, to neglect Dj, S_j_i and Sj+i, then it is neces- 
sary that [I /St] be larger than each of the neglected values at every grid 
point and for all times, i.e., [l/St] > \S_j-i\ + \Sj+i\ Vj G (1,N). This is 
equivalent to require: 

iml ■ Ar m Ar m -, 



IjL-LI j + 1 ~T J — f— 1 j 

In plane geometry, this inequality reduces to: 
AXj AXjAX j+1 



8t <(—±)( 



AXj + AXj+i' 



which must be fulfilled at every finite volume cell. This inequality is 
extremely similar to the CFL condition resulting from von Neumann 
stability analysis. 
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Therefore, the narrow range of stability of explicit methods is a consequence 
of over-simplification of the coefficient matrix. 



3 System of equation - general case 



The set of 2D-hydrodynamical equations in conservative form and in Cartesian 
coordinates may be written in the following vector form: 

^ + L XiXX F + L y , yy G = /, (15) 



where F and G are fluxes of q, and L XjXX , L ym are first and second order 
transport operators that describe advection-diffusion of the vector variables q 

— * 

in x and y directions. / corresponds to the vector of source functions. 

By analogy with Eq.7, when linearization Eq. 15, the following matrix form 

can be obtained: 

[±- + 6(AL X , XX + BL y>yy - H)]Sq = RHS, (16) 



where A[= dF/dq], B[— dG/dq] and H[— df/dq], are evaluated on the former 
(or the new, if Newton method is used to construct the matrices) time levels. 
RHS n = [f — L X:XX F — Ly !yy G] n . 

Adopting a five star staggered grid discretization, it is easy to verify that at 
each grid point Eq. 16 acquires the following block matrix equation: 

^ + _^ j . 1 , t+ ^ J , t+ ^ J+1 , k 

+ £ y 5g jlk -i + D*5q hk + S y S %k+1 = RHSj, k , (17) 



where the underlines (overlines) mark the sub-diagonal (super-diagonal) block 
matrices in the corresponding directions (see Fig. 3). D x ' y are the diagonal 
block matrices resulting from the discretization of the operators L XjXX F, L ytjy G 
and /. 

To outline the directional dependence of the block matrices, we re-write Eq. 
17 in a more compact form: 

+S x 5 qi _ hk +D mod 5q ijk +S x 5q i+ljk = RHS iM (18) 
+SJ5q hk _ 1 . 
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Heat Diffusion 
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Fig. 6. The one-dimensional heat diffusion problem. The profile in (a) corresponds 
to the initial distribution of the temperature. The T-profile in (b) is a steady so- 
lution which was obtained without a source term, i.e., by solving the equation: 
T t = r _2 (r 2 z/T r ) ir (see Sec. 5.1). The profile (c) corresponds to the stationary 
solution obtained with a source term (see Eq. 21). 
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Fig. 7. The one-dimensional heat diffusion problem. The evolution of the 
CFL- number {y St/ Ax 2 , solid line) and the residual (dashed line, RHS of Eq. 7) ver- 
sus number of iteration for different numerical methods. While the explicit method 
stagnates after 130000 time steps, the semi-explicit method converges to the steady 
solution after 20000 time steps. On the other hand, the fully implicit method pro- 
vides the sought solution after 100 iterations only. Although in the latter case the 
defect-correction strategy in combination with several global iterations have been 
employed to assure convergence, the implicit solution procedure is remarkably more 
robust and more efficient than the former methods. 
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where D mod = 5q^ k /5t + D x + D y . The subscripts "j" and "k" denote the grid- 
numbering in the x and y directions, respectively (see Fig. 3). Eq. 18 gives rise 
to three different types of solution procedures: 

(1) Classical explicit methods are obviously very special cases that are recov- 
ered when the sub- and super- diagonal block matrices together with D x 
and D y are neglected. The only matrix to be retained here is (1/St) x (the 
identity matrix), i.e., the first term on the LHS of Eq. 16. This yields the 
vector equation: 



(2) Semi-explicit methods are recovered when neglecting the sub- and super- 
diagonal block matrices only, but retaining the block diagonal matrices. 
In this case the matrix equation reads: 



We note that inverting D mod is a straightforward procedure, which can 
be maintained analytically or numerically. 
(3) A fully implicit solution procedure requires retaining all the block ma- 
trices on the LHS of Eq. 18. This yields a global matrix that is highly 
sparse (a/Fig. 1). In this case, the "Approximate Factorization Method" 
(-AFM: Beam & Warming 1978) and the "Line Gauss-Seidel Relaxation 
Method" (-LGS: MacCormack 1985) are considered to be efficient solvers 
for such a set of HD-equations in multi-dimensions. 

3. 1 Residual smoothing and accelerating convergence 

Let [a,b] be the interval on which Eq. 2 should to be solved. We may divide 
[a,b] into N equally spaced finite volume cells: Ax\ = (b — a)/N, % — 1, N. To 
follow the time-evolution of q using a classical explicit method, the time step 
size must fulfill the CFL-condition, which requires St to be smaller than the 
critical value: 8t™ If [a,b] is divided into N highly stretched finite volume cells, 
for example Ax\ < Ax2--- < Axn, then the CFL-condition restricts the time 
step size to be even smaller than 5t^ n = mmi{(Axi/(V + Vg)i}, which is much 
smaller than St". Thus, applying a conditionally stable method to model flows, 
while using a highly non-uniform distributed mesh, has the disadvantage that 
the time evolution of the variables in the whole domain are artificially and 
severely affected by the flow behaviour on the finest cell. 
Moreover, time-advancing the variables may stagnate if the flow is strongly or 
nearly incompressible. In this case, V$ >> V, which implies that the time step 
size allowed by the CFL-condition approaches zero. 



(19) 



(20) 
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However, we may still associate a time step size with each grid point, e.g., 
5t- m = Axi/(V + Vg), and follow the time evolution of variable q in each 
cell independently. Interactions between variables enter the solution proce- 
dure through the evaluation of the spatial operators on the former time level. 
This method, which is occasionally called the "Residual Smoothing Method" 
proved to be efficient at providing quasi- stationary solutions within a reason- 
able number of iterations, when compared to normal explicit methods (Fig. 
10, 13, also see Swanson & Turkel 1997; Enander 1997). 

The main disadvantage of this method is its inability to provide physically 
meaningful time scales for features that possess quasi- stationary behaviour. 
Here we suggest to use the obtained quasi- stationary solutions as initial con- 
figuration and re-start the calculations using a uniform and physically relevant 
time step size. 



4 Other similar approaches 

The Dufort-Frankel scheme is most suited for modelling diffusion-dominated 
flows and provides solutions of second order spatial and temporal accuracies 
(see Fletcher 1991 for further details). The scheme is unconditionally stable, 
as it relies on adding the positive coefficient resulting from finite difference 
discretization of the diffusion operator to the diagonal elements, thereby en- 
hancing the diagonal dominance of the matrix coefficient. 
Despite this similarity, the scheme differs still from the explicit-implicit meth- 
ods presented here at least in two aspects: 

(1) The DuFort-Frankel scheme is generically second order accurate in time. 
The scheme fails to be consistent when a large time steps are used. Thus, 
the method is not suited for searching steady state solutions. 

(2) The scheme is not suited for modelling advection-dominated flows. 

Another solution procedure that might look similar to the semi-explicit method 
presented here is the Jacobi iteration method. Inspection of Eq. 20, however, 
shows that neither the lower nor the upper diagonal matrices are considered in 
the present solution scheme. Furthermore, the consistency of our solution pro- 
cedure with the original set of equations is guaranteed through adopting the 
defect-correction strategy. Unlike Jacobi method which may diverge if large 
time steps are used, the test calculations presented here show that the semi- 
explicit method converges even for relatively large CFL numbers. 
Recently, Yabe et al. (2001) have presented the Constrained Interpolation 
Profile (CIP) method, which has been successfully applied to model solid, liq- 
uid, gas and plasmas. The method is a kind of semi-Lagrangian scheme which 
can re-produce strongly time-dependent solutions with large CFL-numbers. 
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Fig. 8. The wave propagation problem. The initial density profile is the sinus wave 
(a), which moves right- wards with velocity U = 1. The computed profiles "b" to "g" 
have been obtained after time = 1.966, using 200 finite volume cells. The adopted 
advection scheme is of third order spatial accuracy and second order accurate in 
time. The profiles "b, d, f" have been obtained using the fully implicit solution 
procedure with time step sizes that correspond to CFL=1, 2.5 and 5, respectively. 
The profiles "c, e, g" have been obtained using the semi-explicit approach with 
CFL=1, 2.5 and 5, respectively. 
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Fig. 9. Free-fall of gas onto a non-rotating black hole. Top: 50 density contours 
of a freely-falling gas around one solar-mass black hole. Bottom: the numeri- 
cally-obtained density distribution has a power law that coincides precisely with 
the theoretical solution r~ 3 / 2 . In this calculation, inflow (outflow) at the outer (in- 
ner) boundary have been imposed. 200 strongly-stretched finite volume cells in the 
radial direction and 60 in horizontal direction have been used. 



Therefore, the method is promising and testing its robustness and capability 
at modelling astrophysical phenomena is extremely useful. 
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Fig. 10. Plasma in Free-fall onto a black hole. Similar to Fig. 7, the evolution of 
the CFL-number and the residual versus number of iteration are shown. The solu- 
tion methods used here are: normal explicit (top/left), semi-explicit (middle/left), 
semi-explicit in combination with the residual smoothing strategy (bottom/left), 
semi-explicit using moderate CFL- numbers (top/right), semi-explicit method in 
which the time step size is taken to be a function of the maximum residual (mid- 
dle/right), and finally the fully implicit method (bottom/right). Obviously, the dif- 
ferent forms of the semi-explicit method used here are stable and converges to the 
stationary solution, though at different rates. 
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5 Test calculations 



In the following, we present the results of several test calculations aimed to 
examine the stability, robustness and the capability of the semi-explicit solver 
to treat strongly time-dependent or steady inviscid flows governed by strong 
shocks. 



5.1 The diffusion equation 



Diffusion of heat is generally expressed by a second order differential opera- 
tor. Taking into account this operator in the von Neumann stability analysis, 
it can be shown that it imposes a further limitation on the time step size, 
hence narrowing the stability regime of explicit methods even further. Since 
this limitation may easily lead to complete stagnation of these methods, espe- 
cially when steady state solutions are sought, several modification strategies 
have been suggested. A popular strategy is to treat the diffusion operators 
implicitly. In the absence of advection and other source terms, however, this 
modification is significant, as the scheme is then fully implicit, and therefore 
have the usual drawbacks of implicit methods. 

The semi-explicit method presented here is almost as efficient as normal ex- 
plicit methods. The main difference between these two approaches is that semi- 
explicit methods require additional programming efforts. Specifically, the pos- 
itive entries resulting from finite differencing of the operators involved should 
be carefully selected and added to the diagonal elements of the coefficient 
matrix. 

The purpose here is to apply the semi-explicit method to the heat diffusion 
problem and compare its convergence behaviour with those of explicit and 
fully implicit methods. 

Similar to the other test calculations presented in this paper, the solvers here 
are constructed using routines from the solver package "IRMHD" (Hujeirat 
& Rannacher 2001). The majority of these routines are designed and pro- 
grammed in spherical geometry. Therefore, to use them properly, the equa- 
tions of interest should be reformulated in spherical geometry. Unlike in plane 
geometry, the main drawback of this procedure are the difficulties associated 
with constructing relevant analytical solutions in spherical geometry. On the 
other hand, locating the domain of calculation far from the center, may re- 
duce the difference between these geometries (~ 1/r), and which can be made 
negligibly small. 
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Thus, the one-dimensional diffusion equation in spherical geometry reads: 

— - ——r 2 v— 4 ; (21) 
dt r 2 dr dr 



where T, v denote the variable temperature and the constant diffusion co- 
efficient, respectively. This equation is solved in the radial direction in the 
interval [ri n ,r out ] = [1000,1003], which is sub-divided into 180 finite volume 
cells of equal size. The heat equation is solved using the boundary condi- 
tions: T(r in ) = T(r out ) = 1. Note that the inner and outer boundaries are 
chosen to be very close in order to reduce geometric compression. Further- 
more, to pronounce the difference in the convergence history between the 
explicit, semi-explicit and fully implicit methods, we set v = 10~ 2 and use 
an initial configuration that is very different from the sought steady solution. 
Specifically, the calculations have been initialized using the starting profile: 
T(r,t = 0) = 10e- 10 ( r - r °) 2 , where r Q = (r in + r out )/2 (see Fig.12). 
To first order in (1/r), the steady solution is very close to: T(r,t = oo) = 
-^(r-r in )(r-r in ) + 1. 

In Fig. 6 and 7 the results of three different solution procedures are presented. 
For simplicity, the T-profiles in this figure are plotted versus f = r — 1000. 
Obviously, the semi-explicit method converges to the sought steady solution 
much faster than explicit, but extremely slower than the fully implicit solver. 
Moreover, the method is numerically stable even when extremely large CFL- 
numbers are used. 



5.2 Wave propagation in one dimension 



In an attempt to examine the capability of the method at capturing strongly 
time dependent features, the following one-dimensional wave equation in spher- 
ical geometry has been solved: 

dp 1 d 9 . , 



where p, U denote the density and radial velocity, respectively. The velocity 
assumes the constant value U = 1. To reduce the effect of geometrical com- 
pression, the interval of calculation is taken to be [100 < r < 104]. Similar to 
the previous test problems, different solution procedures have been adopted 
using advection schemes of third order spatial accuracy and second order ac- 
curate in time. 

The initial density distribution, is taken to be a sinus profile that moves right- 
wards with velocity U — 1 ( see Fig. 8, plot a). As in plane geometry, the 



21 



shape of the initial density profile is expected to be preserved as it moves 
right-wards. 

To compare results, a reference solution has been produced using an explicit 
solution procedure. The interval of calculation has been sub-divided into 1000 
finite volume cells, and a time step size that corresponds to CFL = 0.01 has 
been used. The advection scheme here is third order accurate in space and 
second order in time. 

Obviously, the semi-explicit method is capable at capturing the propagation 
of the wave front accurately (see Fig. 8). The amplitude of the wave is hardly 
changed, even when using a time step size that corresponds to a CFL-number 
larger than unity. However, in the latter case, several iterations per time step 
were required to assure convergence. For example, the profile "e" of Fig. 8 
has been obtained using CFL= 2.5 in combination with two global iterations 
per time step, which apparently are not sufficient for convergence. However, 
increasing the number of iteration from 2 to 7, a more accurate solutions has 
been obtained even when the CFL number is increased from 2.5 to 5 (see Fig. 
8 /profile "g"). Consequently, convergence can be maintained if a sufficient 
number of global iterations is performed 4 , though the exact correlation is hard 
to deduce. We note, however, that the LHS and RHS of Equations 7 or 16 
may differ significantly from each other if the flow is strongly time-dependent. 
Here the RHS is calculated using higher order spatial and temporal accuracies, 
whereas the LHS is calculated using first order accuracies. Maintaining com- 
patibility, in this case, depends on the number of global iterations performed, 
which is expected to increase with increasing the CFL-number, but which may 
diverge if the flow contains discontinuities or shocks. 

The profile "f" of Fig. 8 show that even fully implicit methods (without iter- 
ation, hence the LHS and RHS are weak compatibile) may fail to reproduce 
the correct wave profile, if large CFL-numbers are used. 

According to Eq. 6, using 9 < 1 yields a certain average of the q— values on 
the new and old time levels. This averaging-process is useful if the sought 
solutions are strongly time-dependent, such as propagations of shock waves, 
where using large time-steps may cause divergence of the solution track, or 
even numerically-destabilize the scheme. A possible way to capture both time- 
dependent and quasi- stationary solutions is to relate 9 to the time step size, 
which in turns depends on the total residual. An example is the damped 
Crank-Nicolson scheme: 9 — (1 + a St) /2 < 1, where a is a constant of order 
unity (see Hujeirat & Rannacher 2001 and the references therein). In the case 
that stationary solutions are of interest, 9 = 1 should be used. 



4 Using high spatial and temporal accuracies may lead to significant deviations of 
the LHS from the RHS of Eq. 7. This has the consequence that several iterations 
might be still required to reduce the RHS to below a certain critical value. 
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5.3 Free-fall onto a black hole 

A non-rotating gas around a spinless black hole is gravitationally bound, and 
therefore should fall-freely onto the black hole, provided that no other forces 
oppose gravity. The density and velocity profiles of a freely-falling matter 
in the radial direction and far from the event horizon obey the power laws: 
r~ 3 / 2 and r~ 1//2 , respectively. This physical problem is relevant for testing the 
capability of the numerical approach at capturing steady and oscillation-free 
solutions, even when a strongly stretched mesh distribution is used. 
The equations to be solved in this problem are the continuity, the radial and 
horizontal momentum equations, and the internal energy equation. The flow 
is inviscid and adiabatic with the adiabatic index 7 = 5/3. The equations have 
been solved using a first order accurate advection scheme both in space and 
time. In carrying out these calculations, the following conditions/inputs have 
been taken into account. 

• The central object is a one solar-mass non-rotating black hole. 

• The outer boundary is 100 times larger than the the inner radius, i.e., i? Q ut = 
100 x R m , where R m is taken to be the radius of the last stable orbit 5 -Rls 
. To first order in V/c, the flow at this radius can be still treated as non- 
relativistic, though the error can be as large as 30%. 

• Along the outer boundary, the density and temperature of the gas assume 
uniform distributions, and flow across this boundary with the free-fall ve- 
locity. Symmetry boundary conditions along the equator, and asymmetry 
boundary conditions along the axis of rotation have been imposed. Along 
the inner boundary, we have imposed non-reflecting and outflow conditions. 
This means that up-stream conditions are imposed, which forbid informa- 
tion exterior to the boundary to penetrate into the domain of calculations. 
In particular, the actual values of the density, temperature and momentum 
in the ghost zone r are erased and replaced by the corresponding values 
in the last zone, i.e, the zone between i? in and R m + AR. In the case that 
second order viscous operators are considered, care has been taken to assure 
that their first order derivatives across R in are vanished. 

The above set of equations are solved in the first quadrant [1 < r < 100] x [0 < 
9 < 7r/2], where 200 strongly stretched finite volume cells in the radial direc- 
tion and 60 in the horizontal direction are used. 

Fig. 9 shows the 2D distribution of the density around the hole. The bottom 



5 i?Lg = 3 x i? s = 6 x i? g , where R§ and R g are the Schwarzschild and gravitational 
radii, respectively. 
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plot of Fig. 9 shows the precise agreement of the numerical with the theoreti- 
cal solution. 

As in Fig. 7, we show in Fig. 10 the evolutions of the CFL-number and the 
residual as function of the number of iteration which has been obtained us- 
ing different numerical approaches. The results show that the convergence of 
the explicit and semi-explicit methods are rather slow when a relatively small 
time step size is used. This implies that the amplitude-limited oscillations are 
strongly time- dependent that may result from geometric compression. Indeed, 
these perturbations disappear, when relatively large time-step sizes are used 
(see Fig. 10, middle/right). 

In addition, the semi-explicit solver has been tested in combination with the 
residual smoothing strategy. As expected, this approach accelerates conver- 
gence considerably ( Fig. 10: compare the plots bottom/left with the top/right). 

In most of the cases considered here, the time-step size is set to increase in a 
well-prescribed manner and independent of the residual. However, determining 
the size of the time step from the residual directly did not provide satisfactory 
convergence histories (Fig. 10, middle/right). 

The results obtained here indicate that the semi-explicit method is stable when 
compared to normal explicit methods. Furthermore, the semi-explicit method 
can be applied to search for stationary solutions using large time steps, or 
equivalently, CFL-numbers that are significantly larger than unity (Fig. 10, 
middle/left). 



5.4 Shock formation around black holes 

Similar to the forward facing step in CFD, a cold and dense disk has been 
placed in the innermost equatorial region: [1 < r < 10] x [—0.3 < 9 < 0.3] 1 . 
We use the same parameters, initial and boundary condition use in the previ- 
ous flow problem. A vanishing in- and out-flow conditions have been imposed 
at the boundaries of this disk. The gas surrounding the disk is taken to be 
inviscid, thin, hot and non-rotating. Thus, the flow configuration is similar to 
the forward facing step problem usually used for test calculations in CFD. The 
disk here serves as a barrier that forbids the gas from freely falling onto the 
black hole, and instead, it forms a curved shock front around the cold disk. The 
purpose of this test is mainly to examine the capability of the method at cap- 
turing steady solution governed by strong shocks. In solving the HD-equations 
(see Sec. 4.1), an advection scheme of third order spatial accuracy and first 
order accurate in time has been used. The domain of calculation is sub-divided 
into 200 strongly-stretched finite volume cells in the radial direction and 60 in 



1 In applying spherical geometry, the transformation 6 = tt/2 — 6 has been used. 
This is useful if further transformation into cylindrical geometry is planned. 
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Fig. 11. Free-fall of gas onto a black hole surrounded by a cold disk. Top: the density 
distribution (red: large values, blue: low values, green: intermediate values). Middle: 
the temperature distribution (red: large values, blue: low values, gray: intermediate 
values). The curved shock front, where the temperature attains maxima is obvious. 
Bottom: the distribution of velocity field. 
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horizontal direction. In Fig. 11 the configuration of the steady distributions 
of the density, temperature and the velocity field are shown. Similar to the 
calculation in the previous sections, the results indicate that the semi-explicit 
method is stable and converges to the sought steady solution even when a 
CFL-number of order 200 is used (see Fig. 12). However, the method con- 
verges relatively slowly compared to the implicit operator splitting approach, 
where steady solutions have been obtained after one thousand iterations only. 
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Fig. 12. Free-fall of gas onto a black hole surrounded by a disk. As in Fig. 7, 
the evolution of the CFL-number and the residual versus number of iteration are 
shown. The semi-explicit method is apparently stable and converges to the sought 
stationary solution, even when large time step sizes that corresponds to CFL-number 
of 220 are used. On the other hand, 1000 iterations were sufficient to obtain the same 
solution using the implicit operator splitting approach (see Hujeirat & Rannacher 
2001). 
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5. 5 Weakly incompressible flows between two concentric spheres 

To verify the capability of the method at modelling weakly incompressible 
flows, Taylor flows between two concentric and rotating spheres has been 
tested (Fig. 14). This is an ideal test problem, as the flow here is viscous 
(second order diffusive operators are included), closed boundaries (without 
external perturbations) and the corresponding equations accept strict time- 
independent solutions. 

Fig. 13 shows the time-development of the CFL-number and the total resid- 
ual for 5-different solution procedures for searching steady state configurations 
for Taylor flows between two concentric spheres. Using spherical geometry, the 
set of the 2D axi-symmetric Navier-Stokes equations are solved. The set con- 
sists of the three momentum equations, the continuity and the internal energy 
equations. The flow is assumed to be adiabatic. 

As boundary conditions, the inner sphere is set to rotate with Q = 5, whereas 
the outer sphere has Q = 0. The density, temperature, radial and angular 
velocities are set to be symmetric along the equator and along the polar axis, 
except the angular velocity which assumes anti-symmetric conditions along the 
axis of rotation. For the horizontal velocity, anti-symmetric conditions are im- 
posed both along the equator and along the axis of rotation. On the inner and 
outer radius of the spheres, all velocity-components are set to vanish. Initially, 
the flow between the two spheres has zero poloidal and toroidal velocities, so 
that the rotational energy is injected into the flow through viscous interaction 
with the inner boundary. In these test calculations, the viscosity coefficient 
7] = 10~ 2 and the switch on/off parameter 0=1 have been used. In the ex- 
plicit case, the equations are solved according to Eq. 19. For the semi-explicit 
procedure, we solve the HD-equations using the block matrix formulation as 
described in Equation 20. The implicit operator splitting approach is based in 
solving each of the HD-equations implicitly. Here, the LGS method has been 
used in the inversion procedure of each equation (Hujeirat & Rannacher 2001). 
Unfortunately, while this method has been proven to be robust for modeling 
compressible flows with open boundaries, it fails to achieve large CFL-numbers 
in weakly incompressible flows (Fig. 13). This indicates that pressure gradi- 
ents in weakly incompressible flows do not admit splitting, and therefore they 
should be included in the solution procedure simultaneously at the new time 
level. 

Using the semi-explicit solution procedure, the CFL-numbers obtained in the 
present modelling of Taylor flows are, indeed, larger than unity (middle, Fig. 
13), but they are not impressively large as we have predicted theoretically. We 
may attribute this inconsistency to three different effects: 1) The flow consid- 
ered here is weakly incompressible. This means that the acoustic perturbations 
have the largest propagational speeds, which require that all pressure effects 
should be included in the solution procedure simultaneously on the new time 
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Fig. 13. The development of the CFL-number (left axis, solid line) and the total 
residual (right axis, dashed line) versus covered-time in normalized units of five 
different numerical methods (from to top to bottom: normal explicit, residual 
smoothing, semi-explicit, implicit operator splitting (Hujeirat & Rannacher 2001) 
and the fully implicit method). While the effective time covered in each run of 
these different methods is similar, the actual number of iteration is substantially 
different. The numerical problem here is to search stationary solutions for Taylor 
flow between two concentric spheres. The inner sphere has a radius rj n = 1 and 
rotates with angular velocity fij n = 5 , whereas the outer sphere is non-rotating and 
its radius is taken to be r ou t = 1.3. 6 = 1 and r? = 10~ 2 have been used in these 
calculations. The initial density and temperature are taken to be p(r, 6, t = 0) = 1 , 
and T(r, 6,t = 0) = 10 1 , respectively. The computational domain is [1, 1.3] x [0, tt/2] 
and consists of 30 x 50 non-uniformly distributed tensor-product mesh. 
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Fig. 14. Steady state solutions of the Taylor flow between two concentric spheres. 
Here the velocity field and the angular velocity (large-to-low values correspond 
to blue-to-red colors) are shown. The capability of the methods to capture the 
formation of rotationally-driven multiple vortices near-equatorial region is obvious. 

level. Indeed, Fig. 15 shows that solving the angular momentum only may re- 
lax these limitations and that extremely large CFl-numbers can be obtained. 

2) The conditions imposed on the boundaries are non-absorbing, and do not 
permit advection of errors into regions exterior to the domain of calculation. 

3) The method requires probably additional improvements in order to achieve 
large CFL-numbers. This could be done, for example, within the context of 
the "defect-correction' iterative procedure, in which the block diagonal matrix 
Anod is employed as a pre-conditioner. 

In the final case (bottom, Fig. 13), the whole set of HD-equations is solved 
taking into account all pressure terms in a fully implicit manner. Here we 
use the AFM for solving the general block matrix-equation which is locally 
described by Eq. 18. 

For controlling the time step size in these calculations, we have adopted the 
following description: St = a e/ max(RHSj^), where e = min(AXj, AY k ), and 
where «o is a constant of order unity. The maximum and minimum functions 
here are set to run over the whole number of grid points. 
Worth noting is the difference between the terminal values of the residual in 
the semi-explicit and the fully implicit cases. After a certain number of itera- 
tion, the residual in the former case does not decrease, and instead it shows 
strong, but still limited variations with increasing the number of iteration. 
Although a CFL-number of order 10 is achieved in this test calculations, the 
origin and magnitudes of these variations are not completely clear. However, 
taking into account that these strong variations do not show up in the fully 
implicit case, we may conclude that their origin is connected to the incomplete 
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Fig. 15. The semi-explicit method applied to a scalar problem. As in Fig. 7, this Fig- 
ure shows the development of the CFL-number (left axis, solid line) and the residual 
(right axis, dashed line) of the angular momentum equation in two dimensions ver- 
sus the number of iteration. As initial conditions we use the steady distributions of 
the physical variables that have been obtained from the simulations of the Taylor 
problem (see Fig. 14). This includes the velocity field, density, temperature and rj. 
For the angular velocity we use Q = as initial condition. 




Table 1 

The three different numerical approaches and their properties 





Explicit 


Semi-Explicit 


Implicit 


Stability: 


stable within 
CFL < 1 


absolutely stable 


absolutely stable 


Convergence 
speed 1 : 


limited by 
CFL < 1 


faster than explicit, 
slower than implicit 


unlimited by the 
CFL condition 


Efficiency 2 : 


highly efficient 


efficient 


inefficient 


Robustness: 


very limited 


robust, but can be made 
highly robust 3 


robust 


Programming 
efforts: 


easy 


easy 


relatively difficult 



inclusion of the source terms in the coefficient matrix. 



6 Summary 



In this paper a strategy for modifying explicit methods and constructing a 
class of semi-explicit methods is presented. We have shown that this modifi- 
cation enhance the robustness and enlarge the range of application of explicit 
methods, and in particular, it enables their use for searching quasi-stationary 
solutions for the Euler and Navier-Stokes equations in multi-dimensions. 
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The new strategy is based on reformulating explicit methods in matrix form, 
i.e, transforming explicit methods from scalar into matrix-vector problems. 
Thus, solving the set of time- dependent HD or MHD equations would require 
solving a matrix equation of the form: Aq = b. Explicit methods in this case 
relies on approximating and replacing the matrix A by the most simple matrix 
in algebra: the identity matrix I. This solution-procedure may converge, if the 
time step size is sufficiently small, or equivalently, if the entries of the matrix 
I/5t are sufficiently large, so that all off-diagonal elements of A can be safely 
neglected. As we have shown in Sec. 2, this yields a condition similar to the 
requirement that CFL < 1. 

Unlike normal explicit methods, in which inclusion of diffusion generally causes 
further limitations on the time step size, diffusion operators in the semi-explicit 
formulation pronounce the diagonal dominance and enhances the stability of 
the inversion procedure, irrespective of the dimensionality of the problem. 
However, the convergence rates of semi-explicit methods remain lower than 
that of their fully implicit counterparts. 

A relevant question which arises here is: how efficient are semi-explicit meth- 
ods compared to explicit, and if they are more favorable than explicit methods 
when modelling time- dependent flows? 

We note that in applying the semi-explicit method to a scalar equation, ap- 
proximately N additional algebraic operations are required for calculating the 
diagonal elements, and roughly N division operations. Thus, applying the 
semi-explicit approach to system of equations of M variables, the additional 
algebraic operations required scale as: 2 x iV x M. If we take into account 
the addition global iterations required per each time step to assure conver- 
gence, we conclude that the semi-explicit methods are, indeed, less efficient 
than their explicit counterparts, provided the flow is strongly time-dependent. 
On the other hand, since the additional algebraic operations increases linearly 
with the total number of variables, the associated computational load is likely 
to be small compared to calculating the RHS of Eq. 7 or 16. Furthermore, 
the algorithm can be optimized further, in such a manner that the relevant 
routines, in which the corresponding additional operations are performed, are 



This applies for diffusion dominated flows, and within the regime of stability of 
the method. The discretization method used is assumed to be consistent /compatible 
with the continuous formulation. Semi-explicit methods converge relatively slowly, 
though not directly affected by the CFL-condition. 

2 A measure for the total computing operations per time step, provided the flow 
is time-dependent. This measure may reverse if the total number of computing 
operation is of interest, and if quasi-stationary solutions are sought. 

3 Reducing the semi-explicit method into explicit when modelling strongly time- 
dependent flows, and enhancing the degree of implicitness when modelling radia- 
tive, MHD and/or weakly/strongly incompressible flows through incorporating ad- 
ditional non-diagonal entries (see Fig. 1). 
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switched off if CFL < 1. In this case, the semi-explicit approach degenerates 
into an explicit method when the underlying flow is strongly time-dependent. 
This optimization procedure is particularly useful when modelling propagation 
of waves and moving shocks. The latter phenomena are non-linear features of 
HD-flows, and the corresponding equations do not admit stationary solutions. 
In this case, the time step size should be chosen sufficiently small in order to 
capture shock profiles accurately, irrespective of whether the method used is 
explicit or implicit. For modelling such flows, explicit methods are superior 
over implicit, provided that the flow is isothermal, adiabatic or polytropic. 
In the case that other physical and chemical processes are concerned, which 
generally operate on much shorter time scales than the dynamical time scale, 
modifying explicit methods into semi-explicit or strongly implicit is unavoid- 
able. 

In Table 1, the main properties of the new semi-explicit method, compared to 
explicit and implicit methods are listed. 

Additionally, we have shown that the residual smoothing approach accelerates 
the convergence of both explicit and semi-explicit methods. 
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